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Abstract 

^ I Based on a novel point of view on 1-dimensional Gaussian quadrature, we present a 

5^ I new approach to the computation of d-dimensional cubature formulae. It is well known 

that the nodes of 1-dimensional Gaussian quadrature can be computed as eigenvalues 
of the so-called Jacobi matrix. The d-dimensional analog is that cubature nodes can 
be obtained from the eigenvalues of certain mutually commuting matrices. These are 
obtained by extending (adding rows and columns to) certain noncommuting matrices 
Ai, . . . , Ad, related to the coordinate operators xi, . . . , Xd, in R*^. We prove a cor- 
respondence between cubature formulae and "commuting extensions" of Ai, . . . ,Ad, 
satisfying a compatibility condition which, in appropriate coordinates, constrains cer- 
tain blocks in the extended matrices to be zero. Thus, the problem of finding cubature 
formulae can be transformed to the problem of computing (and then simultaneously 
diagonalizing) commuting extensions. We give a general discussion of existence and 
of the expected size of commuting extensions and describe our attempts at computing 
them, as well as examples of cubature formulae obtained using the new approach. 



1 Introduction 



One of the most elegant topics in numerical analysis is the theory of Gaussian quadrature 
[1] . Unfortunately this theory is limited to one dimension, and although something is known 
about generalizations to multiple dimensions (see [2] for a survey article and many refer- 
ences), at the moment there are many more questions than answers. The aim of this paper is 
to present a new approach to cubature rules ( "cubature" seems to be the name given to the 
generalization of quadrature to arbitrary dimension). In classical, one-dimensional Gaus- 
sian quadrature, the most widely used method for computing nodes and weights, developed 
about 35 years ago [3], involves solving the eigenproblem for a certain tridiagonal matrix 
(see [4] for a recent "basis independent" discussion of this). As far as we are aware, no ex- 
tension of this to higher dimensions has previously been obtained. Our proposed method for 
computing rf-dimcnsional cubature formulae involves the construction of d matrices (with 
tridiagonal block structure in suitable bases), extending these, in a manner we will explain 
below, to a set of commuting matrices, and then solving the simultaneous eigenproblem for 
these commuting matrices. 

The main novel step in this process, which follows very naturally from a new approach 
we present to one-dimensional Gaussian quadrature, is the need to construct commuting 
extensions of a set of matrices. We say the N x N matrices Ai, A2, . . . , are N x N 
commuting extensions of the n x n matrices Ai, A2, . . . , A^ (here N > n) ii the top left 
n X n block in Aj is A^, for each i — 1, 2, . . . , d, and the matrices ^i, A2, . . . , pairwise 
commute. The idea of commuting extensions is very natural, but we do not find any such 
notion in the linear algebra or numerical linear algebra literature (see, however, some very 
similar ideas in a recent, independent, poster of Cargo and Littlejohn [5]). Since we hope 
the idea of commuting extensions will find other applications, the first few sections in this 
paper explore this subject without reference to cubature rules. Section 2 covers basic theory, 
section 3 discusses a couple of simple algorithms for computing commuting extensions, with 
which we currently have very limited success, and section 4 discusses commuting extensions 
when the matrices A^ take a special form, relevant for the study of cubature rules. 

In section 5 we turn to the theory of cubature rules. Subsection 5.1 contains a novel 
approach to one dimensional Gaussian quadrature, based upon the properties of a certain 
operator, its eigenvalues and eigenf unctions. This serves as the model for all the subsequent 
discussion. In subsection 5.2 we consider the natural extension of this approach to multiple 
dimensions, and prove the central results of the paper, giving an equivalence between odd 
degree, positive weight cubature rules and commuting extensions (satisying a compatibility 
condition that will be explained in the sequel) of a certain set of matrices. Subsection 5.3 
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includes some simple consequences of this relationship. One of the key results in the theory 
of cubature rules, a lower bound on the number of nodes needed for an odd degree cubature 
rule, originating in the work of MoUer [6], follows from a general result in the theory of 
commuting extensions (theorem 2 in section 3). Similarly a simple result on the spectra of 
commuting extensions (theorem 6 in section 3) gives interesting constraints on the nodes in 
positive weight cubature rules, which we believe have hitherto been overlooked even in one 
dimension. 

In section 6 we turn to actual application of the commuting extension approach for 
computing nodes and weights. As in section 3, our achievements here are rather limited, 
nevertheless they validate our approach, and even give a few new cubature formulae. Section 
7 contains a list of open questions. 

We close this introduction by mentioning that one of the oldest continuing applications 
of Gaussian quadrature is in quantum mechanics, where it is used, in the so-called "DVR 
method" for the computation of matrix elements of non-exactly solvable Hamiltonians (see 
[7] for early references and [8] for recent reviews). This paper was born out of an attempt 
to extend the DVR method to higher dimensions; other, recent progress on this subject has 
been made by Dawes and Carrington [9]. Another interesting perspective on DVR can be 
found in the paper [10]. 



2 Commuting Extensions 



In this section we present the basic theory of commuting extensions. 

Definition. We say the N x N matrices Ai, A2, . . ■ , are N x N commuting extensions 
of the n X n matrices Ai, A2, . . . , (here > n) if the top left n xn block in A^ is Ai, for 
each i — 1,2, . . . ,d, and the matrices Ai, A2, . . . , A^ pairwise commute. 

Theorem 1. Any set of matrices admits commuting extensions. 

Proof: We construct explicit commuting extensions of the n x n matrices Ai, A2, . . . , A^. 
Take 



Ai A2 As ... Ad 
Ad Ai A2 ... Ad-i 
Ad-i Ad A ■■■ Ad-2 



\ 



As Aa 



^ A2 A3 A4 ... Ai 
Ai A2 A3 ... Ad 
Ad A A2 ... Ad-i 



A3 Aa Ar, 



\ 



etc. (1) 
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More fully, take to he adnx dn matrix which is a dx d matrix oi nxn blocks, with the 
j, kth block, which we denote {Ai)jk, equal to Ai+^^j (mod d)- Then 

d _^ _^ d 

iAiAii)jj/ — = Aj^k-j (mod d)Ai'+j'-k (mod d) ■ (2) 

k=l k=l 

Replacing the summation index k by k' + i' — i (mod d) the second sum becomes 

d 

5Z ^i'+k'-j (mod d)^i+j'-k' (mod d) (3) 

fc'=l 

which is equal to {AiiAi)jji. Thus 

(iiAO^i' = (4) 

for all i.e. and A^r commute. • 

Theorem 1 establishes the existence of commuting extensions, but it is natural to ask 
what is the smallest possible dimension for commuting extensions of a given set of matrices. 
To this end we have the following result: 

Theorem 2. If N x N commuting extensions of the nxn matrices Ai, A2, . . . , A^ exist, 
then 

Ar>n + lmax,,rank([A„.l,]). (5) 

Proof: Suppose Ai, A2, . . . ,A(i are N x N commuting extensions of the nxn matrices 
Ai,A2,...,Ad. Write 

\ bi ai J 

where the matrices Oj, 6j, ccj have sizes nx (N — n), (N—n) xn, (N—n) x (N—n) respectively. 
The top left nxn block of the equation [Ai, Aj] — gives the requirement 

[Ai, Aj] + Gibj - Qjbi = . (7) 

Since the matrices Oj and bi do not have rank exceeding {N — n), neither do products of the 
form Qibj, and the matrices Uibj — ttjbi can have rank at most 2{N — n). Thus (7) can hold 
only if for each i,j we have 

T8ink{[Ai,Aj\) <2{N -n) , (8) 

and the theorem follows directly. • 

Unfortunately there is a large gap between the lower bound on N from theorem 2 and 
the N in the existence proof of theorem 1. In practice, it seems that the lower bound of 



theorem 2 is rarely attained, and the N of theorem 1 is much too big. As we shall see in 
Section 5, theorem 2 gives rise to a well known lower bound on the number of points needed 
for a cubature formula, and in that context also the bound can rarely be attained. 

In addition to not knowing, in general, any way to rigorously predict the lowest dimension 
for commuting extensions of a given set of matrices, we also currently have no way of 
determining how many distinct families of commuting extensions of a given dimension exist. 
By a family we mean a set of commuting extensions related by conjugation as described in 
the following obvious result: 

Theorem 3. If the matrices Ai, A2, . . . , are N x N commuting extensions of the nxn 
matrices Ai, A2, . . . , A^, then so are the matrices UAiU~^, UA2U~^, . . . , UAJJ'^, where U 
is any matrix of the form 



with U an invertible {N — n) x [N — n) matrix. 

To proceed further, and at least get some idea of the size needed for commuting exten- 
sions, we have to resort to parameter counting. From here on we restrict to the case where 
the matrices Ai and Ai are symmetric, i.e. the case of symmetric commuting extensions of 
a set of symmetric matrices. Note that except when d — 2 the existence construction of 
theorem 1 does not guarantee symmetric commuting extensions. Neither is it clear that the 
lowest dimension commuting extensions of a set of symmetric matrices need necessarilly be 
symmetric. But because the case of symmetric commuting extensions of symmetric matrices 
is relevant for cubature rules, we restrict our attention to this. 

If the matrices A^ are symmetric then we can write 



where is n x {N — n) and ai is {N — n) x {N — n) and symmetric. Thus the number of 
free parameters we have in choosing the extensions of the A^ is 




(9) 




(10) 




(11) 



Let us assume that at least one of the A^, say Ai, has distinct eigenvalues. Then all matrices 
that commute with Ai also commute amongst themselves, and we just need to check that 
[Ai, Ai\ = for i = 2, . . . , d. Since the commutator of symmetric matrices is automatically 
antisymmetric, we have 
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equations to satisfy. We cannot, however, directly compare the number of parameters from 

(11) with the number of equations from (12), as from theorem 3 we learn that (except when 
N — n + 1) commuting extensions exist in families. For symmetric commuting extensions 
the matrices U (and thus U) in theorem 3 are restricted to be orthogonal. So symmetric 
commuting extensions occur in families with ^{N — n){N — n — 1) parameters, and the 
number of parameters in choosing extensions should exceed the number of equations from 

(12) by at least this amount. Thus we need 

^d{N -n){N + n+l)> ]-N{N - l){d - 1) + - n){N - n - 1) . (13) 

A little rearranging of this inequality gives the condition 

^n{n-\){d-\) d-1 d^-1 d{d''-l) fl\ 

N -n> ^ ^, ,, — - = — h + o - . (14) 

- 2(n + d) 2 2 2n \nj ^ ' 

If N satisfies this condition we expect to find N x N symmetric commuting extensions. For 
comparison, in the explicit commuting extensions of theorem 1 (which, however, were not 
symmetric) we had N — n = {d — l)n; we can clearly expect to do much better than this. 

On occasions it seems that parameter counting can be misleading. As an example 
consider the case of c? = 2. From theorem 2 we have N — n > |rank([Ai, A2]) (note that 
since the commutator of 2 symmetric matrices is antisymmetric, its rank is always even). 
The parameter counting argument tells us that we should expect 

^_ n(n-l)_ 

- 2(n + 2) ^ ^ 

Assuming [Ai , A2] of maximal rank we have 

irank([A,^]) = ( " . (16) 



So by theorem 2 



, „ n even , , 

N-n>{ . (17 

^ nodd 



We see that when d — 2 and [^1,^2] is of maximal rank the inequality from parameter 
counting is actually weaker than the rigorous one from theorem 2. Thus in this case the 
parameter counting argument is certainly fiawed. But this seems to be rather exceptional; in 
general it appears that when it is consistent with the lower bound of theorem 2, parameter 
counting gives a better idea of the size we should expect for commuting extensions. In 
particular we give the following example where the lower bound of theorem 2 cannot be 
attained, at least with symmetric extensions: 
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Theorem 4. For n > 5 there exist symmetric nxn matrices Ai, A2 with rank ([^i, A2]) = 2 
and no (n + 1) x (n + 1) symmetric commuting extensions. 

The proof proceeds through the following lemma which will also be useful later: 

Lemma 1. Let Ai, A2 be a pair of symmetric nxn matrices with rank {[A-i, A2]) — 2. Let 
{v, w} be a basis of Im {[Ai, ^2])- If ^i- ^2 have (n + 1) x (n + 1) symmetric commuting 
extensions then the vectors v, w, Aiv, Aiw, A2V, A2W are linearly dependent. 

Proof (of lemma 1). (n+ 1) x (n+l) symmetric commuting extensions of Ai, A2 must take 
the form 



Finding (n + 1) x (n + l) symmetric commuting extensions of Ai, A2 is equivalent to finding 
vectors a, 6 and scalars a,P satisfying (19)-(20). 

Since rank {[Ai, A2]) = 2, the vectors a, b cannot be linearly dependent, for if they were 
we would have ab^ — baF — 0, giving a contradiction with (19). Thus we can find a vector 
orthogonal to a but not to b. Applying (19) to this we deduce that a is in Im([Ai, A2]). Like- 
wise, applying (19) to a vector orthogonal to b but not to a, we see that b is in Im([y4i, A2]). 
Choose {v, w} to be a basis of Im ([Ai, A2]). From the previous remarks it follows that we 
can write 




(18) 



where a, b are n-dimensional column vectors and a, (5 are scalars. The requirement [Ai, /I2] 
translates into the equations 



[Ai,A2] + ab^ -bJ = 0, 
Alb + (5a - A2a - ab = 0. 



(19) 
(20) 



a — \v -\- jJLW 



b = 



vv + pw 



(21) 

(22) 



where A, p are constants with \p — p,v ^ 0. Substituting in (20) we have 



(/5A — av)v + {fi\x — ap)w + uAiv + pAiw — XA2V — 11A2W = . 



(23) 



Since Xp—jii/ 7^ we see at once that the 6 vectors v, w, Aiv, Aiw, A2V, A2W must be linearly 
dependent. • 

Proof (of Theorem 4). Let 



where the are distinct, the are arbitrary, and the Wa are entries of two arbitrary 
hnearly independent n-dimensional vectors. Then [^i, ^2] = wv'^ —vw^ , so rank {[Ai, A2]) = 

2, and {v, w} is a basis of Im {[Ai, A2]). By the lemma there can only exist (n + l) x {n + 1) 
symmetric commuting extensions of Ai,A2 if the 6 vectors v,w, Aiv, Aiw, A2V, A2W are 
linearly dependent. There is no evident reason why they should be linearly dependent if 
n > 5, but to show concretely that they usually are not, we considered the case n — 6 where 
the Xa take the values 1,2,3,4,5,6, and then constructed, using Maple, the determinant 
of the 6x6 matrix with columns v,w, Aiv, Aiw, A2V, A2W. The values of the Va,Wa,iJ>a 
were not fixed. The resulting determinant, which is too long to reproduce here, is simply a 
polynomial expression in the eighteen variables Va,Wa,fJ.a, and is not indentically zero. The 
case n > 6 follows trivially from the case n = 6. • 

Note on index conventions: In discussion of commuting extensions we start with d 
matrices of size n x n which we extend to size N x N. For clarity, in most of this paper we 
adhere to the following index conventions: 
Indices i,j, k etc. run from 1 to d. 
Indices a, b, c etc. run from 1 to n 
Indices a, /3, 7 etc. run from 1 to N. 

The results up to here all concern the existence and size of commuting extensions. For 
the purpose of finding commuting extensions, to be discussed in the next section, we will 
use the following: 

Theorem 5. The n x n symmetric matrices Ai, A2, . . . , A^ admit N x N symmetric com- 
muting extensions if and only if there exist N x N diagonal matrices Ai, A2, . . . , A(i and an 
n X N matrix Q with orthonormal rows such that 

A = QKQ^ . (25) 

Proof. From the extensions to Ai,Q: If we can find NxN symmetric commuting extensions 
^1, A2, ■ ■ ■ , Aa, then we can find diagonal matrices Ai, A2, . . . , A^; and an N x N orthogonal 
matrix Q such that 

Ai = QAiQ'^ . (26) 

The matrix Q in the theorem is comprised of just the first n rows of Q. 

From Aj, Q to the extensions: A matrix Q as described in the theorem can always be 
extended, by the addition oi N — n orthonormal rows, to a,n N x N orthogonal matrix Q. 
(In fact this can be done in many ways, corresponding to the freedom described in Theorem 

3. ) Once such a Q has been constructed the matrices = QA^Q^ are NxN symmetric 
commuting extensions of the A^. • 
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Note: The matrix Q in theorem 5 satisfies QQ^ — Inxn- 

It is of interest to understand how the spectra of commuting extensions (i.e. the entries 
of the matrices Aj in Theorem 5) are related to the spectra of the original matrices Aj. 
The following is a first result in this direction. It is a simple consequence of the Sturmian 

separation theorem [11], but we offer a direct proof too. 

Theorem 6. Let A be an iV x symmetric extension of the n x n matrix A. Then the 
smallest eigenvalue of A is less than or equal to the smallest eigenvalue of A, and the largest 
eigenvalue of A is greater than or equal to the largest eigenvalue of A. 

Proof. Diagonalizing A and A we have 

A - QAQ'^ , A = UDU^ , (27) 

where Q and U are N x N and n x n orthogonal matrices, respectively, and A and D are 
N X N and n x n diagonal matrices containing the eigenvalues of A and A, respectively. 
Restricting the first equation to its upper n x n block we have also 

A = QAg^ , (28) 

where Q is an n x A^ matrix composed of the first n rows of Q, which are orthonormal. 
Comparing the two formulae for A and using the orthogonality of U we have 

D = iU^Q)AiU^Qf . (29) 

The matrix U^Q, like Q, has orthonormal rows; in particular we have 

N 

E(^^Q)L = 1, a = l,...,n. (30) 

a=l 

The diagonal entries of (29) read 

N 

Da=J2iU^Q)LK, a=l,...,n. (31) 

a=l 

Thus each eigenvalue Da of A is a generahzed average of the eigenvalues {Aq,} of A with 
respect to a set of positive weights {{U^Q)la}- follows at once that it is not possible that 
any of the Da be either smaller than, or greater than, all of the Aq. • 

3 Computing Symmetric Commuting Extensions 

The most obvious approach to computing commuting extensions is simply to treat the 
unknown entries in the extended matrices Ai as variables, and to consider the conditions 



9 



= as equations in these variables. In the generic case (generically we should 
expect the Ai have distinct eigenvalues) it will be sufficient to look at the equations just for 
one particular value of i. li N — n > 1, then by theorem 3 we expect continuous families 
of extensions, this freedom can be exploited to fix some of the variables. The system of 
equations we obtain will be quadratic in the unknown variables, and can be tackled by 
standard methods for solving systems of equations. We have done some initial experiments 
with this approach in the case d — 2, attempting to solve the system of quadratic equations 
1) by integrating the gradient fiow v' — —V\\[Ai{v), A2{v)]\\'^ (here v denotes the variables 
added to form the commuting extensions), and 2) using Newton's method. The results are 
very varied; for some pairs of moderate-sized matrices there is reasonable convergence, but 
in other cases there are signs of extreme ill-coniditioning (very low gradients in the case of 
gradient flow, almost singular Jacobian in Newton's method). Some of the cubature related 
results we will present in section 6 were obtained by the integration of the gradient flow 
with Frobenius norm of the commutator; some technical details of these calculations can be 
found in [12]. 

In this section we focus on a different approach to computing commuting extensions, 
based on theorem 5 from section 2, and thus restricted to the case of symmetric extensions 
for symmetric matrices. We attempt to construct the matrices Aj and Q introduced in the 
theorem by minimization of 



= E {{A-QKQ^)a,y . (32) 

^ 1=1 a,b=l 

Consider flrst the variation of this with respect to the entries of Aj (which is diagonal) . It 
is straightforward to check that 

■^{Q^{A-QA,Q^)QU { ^ , (33) 



diA,] 

implying that the entries of Aj must satisfy 




^ \ i^l d 

i:{Q^Q)UKh,-^{Q^A,QU ■ (34) 

/3=i [ a = l,...,N 

We see that if we can construct Q, we can, using the last formula, also construct Aj, at least 
assuming invertibility of the matrix with entries {Q^QYafi- 

To build Q, or more precisely Q, the N x N orthogonal matrix that is an extension of 
the nx N matrix Q by addition of {N — n) more rows, we use a sequence of Jacobi rotations 
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(not using the rotations that just mix the last {N — n) rows). More fully: assuming we have 
an initial guess Qo for the matrix Q, we construct a new guess in the form Q — R{9)Qq 
where R{9) is a Jacobi rotation of 2 specified rows of Qq through angle 9. 9 is chosen to 
minimize S, where in S we use equation (34) to determine the Aj. The minimization is done 
numerically, there does not seem to be an explicit formula for 9 available. 

Note that this algorithm will actually solve a more general problem than that of finding 
commuting extensions. Suppose wc use a value of N for which there are no commuting 
extensions. We can still compute a minimum of S, which means we will find Q and Aj 
such that Ai ^ QKiQ^ . Extending Q to an iV x iV orthogonal matrix Q we will have 
commuting matrices = QAiQ^ which are approximately extensions of the A^, i.e. we 
will be computing commuting approximate extensions. In the case N = n we will be 
computing commuting approximations to the matrices Ai without increasing the dimension. 
The question of the existence and computation of commuting matrices which approximate 
a given set of matrices, with commutators of small norm, has a substantial mathematical 
history, having first been asked, apparently, in [13]. For some recent references see [14]. 
For the case of commuting approximations {N — n), the numerical approach we have just 
outlined was given in [15], though unlike in our more general case, it seems that in this case 
there is an explicit formula for the choice of rotation angle 9 at each stage. This method is 
also advocated for numerical simultaneous diagonalization of commuting matrices. There 
are applications of the algorithm in statistics and in signal processing [16], and it is also 
employed in [9]. 

It remains to report on some numerical experiments. In all cases we implemented the 
algorithm starting with various random choices of Q and applying sweeps of Jacobi rotations 
of all possible pairs of rows. Wc attempted to compute symmetric commuting extensions 
(and approximate extensions) for 4 specific pairs of symmetric 6x6 matrices; the findings 
are supported by many other numerical experiments as well. The specific matrices used can 
be found on the internet at 

http : //www. math. biu. ac . il/~schif f /commext .html . 
The pairs of matrices chosen covered the following 4 cases: 

1. rank([Ai, A2]) — 2, 7 x 7 symmetric commuting extensions known to exist. 

2. rank([74i, ^2]) = 2, 7 x 7 symmetric commuting extensions known not to exist (see 
lemma 1, section 2). 

3. rank([74i, A2]) = 4, 8 x 8 symmetric commuting extensions known to exist, this being 
the lowest dimension allowed by theorem 2, section 2. (In fact presumably many 
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commuting extensions exist, as even after accounting for the invariance of theorem 3, 
section 2, there are more parameters than equations). 

4. rank([Ai, A2]) = 6, 9 x 9 symmetric commuting extensions known to exist, this being 
the lowest dimension allowed by theorem 2, section 2. (Once again presumably many 
commuting extensions exist.) 

In each case we not only ran the algorithm with N large enough that we expected to 
find commuting extensions, but also with all other > 6 smaller than this (thus computing 
commuting approximations for = 6 and commuting approximate extensions for N > Q 
but too small for commuting extensions). Our observations can be summarized as follows: 



1. For the actual computation of commuting extensions {N = 7 in case 1, = 8 in cases 
2 and 3, = 9 in case 4) the positive result is that in each case the algorithm was 
observed to converge. After initial stabilization the value of S was observed to drop 
off according to 

hiS^a-bk, (35) 

where k is the sweep number, and a and 6 > are (case-dependent) constants. This 
behavior was observed over many orders of magnitude. The problematic result is that 
except in case 2 the value of h was found to be small, to the extent that obtaining 
even single precision accuracy required many thousands of sweeps. This behavior 
is not particularly surprising given that wc are in practice doing a multidimensional 
optimization by sequential one-dimensional optimizations, with a fixed choice for the 
search directions. Remarkably in case 2 a large value of h was observed (and this 
persisted for other choices of matrices with rank 2 commutator but for which we were 
searching for extensions with N = 8). 

2. For the computation of commuting approximations {N — n), in all cases the algorithm 
converged reasonably quickly, as reported in the literature [15], though there were 
noticeable differences in the rate of convergence, with case 1 being substantially slower 
than all the other cases. No problems with local minima were observed, different initial 
choices of Q produced the same minimum value of S. Testing of this was hmited, 
though. 

3. In the search for commuting approximate extensions {N = 7 in cases 2 and 3, N — 7,8 
in case 4) a variety of behaviors were observed. For A^ = 7 in case 2 the search was, by 
far, the slowest performed, but it did ultimately converge. In case 4 the N — 8 search 
was noticeably slower than the N — 7 search, and for A^ = 8 we noticed convergence 
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to different values of S, indicating the presence of local (approximate) minima. In all 
cases the values of S achieved were at most a few orders of magnitude smaller than 
the values obtained when N — n. This gives hope that in the general case, when we 
do not know what the smallest dimension for commuting extensions is, we might be 
able to detect it by minimizing S for different N, but a better minimization algorithm 
than the current one will certainly be necessary. 

It is clear from our results that a lot more work is necessary on the topic of computing 
commuting extensions. 



4 A Special Case Of Commuting Extensions 



We turn to consideration of a special case of commuting extensions that turns out to be 
relevant for cubature formula. It is characterized by 2 conditions: First, the symmetric nxn 
matrices Ai for which we wish to find commuting extensions have tridiagonal block form 



A, 



\ 



an 


an 


. 










T 

ail 




«i2 ■ 













T 

ai2 


ttiS ■ 
















. 




-1) 


ai{r~l) 








. 


T 


-1) 





\ 



(36) 



Here an, Oii2, . . . ,air are symmetric square matrices of sizes rii x rii, 712 x 712, . . . , x 
respectively, where rii + n2 + ■ ■ ■ + rir — n. The matrices Oji, cii2, ■ ■ ■ , ai(r-i) are of size 
ni X ^2,^2 X n^-i ■ ■ ■ ,nr-i x rir respectively. The second condition we impose is that the 
commutator matrices [Aj, Aj] all vanish except for a single block in the bottom right hand 
corner, of size x n^. 

Let us seek symmetric commuting extensions with the matrices Ai, of size N x N, also 
taking tridiagonal block form, that is 

/ 



an Oil 


. 















(Ii2 ■ 
















OiiS ■ 
















. 


ai(j-- 


-1) 


ai(r-l) 








. 


T 


1) 


aif 







. 







T 

at 





(37) 
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where the new blocks ctj are of size {N — n) x {N — n) and are symmetric, and the are 
of size rir X {N — n). 

The questions we wish to ask are (1) what are the equations that the new blocks aj 

have to satisfy? and (2) how large need N be for us to have a hope that such extensions 
exist? As in section 2, we assume that Ai has distinct eigenvalues, so we need only check 
that Ai commutes with the d — 1 matrices A2, ■ ■ ■ ,Aii and this guarantees that all the Ai 
mutually commute. A brief calculation, using the fact that the commutators [74i,74j] are 
zero except for a single block, gives the following conditions: 



ai(r_i)aj — ai(r-i)ai 

rj-i rp rj-i rj-i 

Cll(r_l)ai(r-1) ~ '^i{r-l)'^i{r-l) + O-XrO-ir — airO-ir + OlOj — OiOi 

a^ai — ajai + aictj — a^ai 










2,...,d 



(38) 

Of these four equations for each i, the first is of size n^-i x {N — n), the second is of 
size rir x rir and antisymmetric, the third is of size Ur x {N — n) and the fourth is of size 
[N — n) X [N — n) and antisymmetric. Thus there is a total of 



{d - 1) [nr-\{N - n) + ^^^(n^ - 1) + rir{N - n) + ]^{N - n){N - n - 1)^ 
equations to be satisfied. The number of variables available in the and ctj is 



d { nr{N -n) + ]^{N -n){N -n+l\ 



The system of equations (38), has an invariance 

ai QiQ , 



ai 



(39) 



(40) 



(41) 



where g is an {N — n) x [N — n) orthogonal matrix. Thus it is not sufficient that the number 
of variables simply exceed the number of equations to be solved, it must exceed the number 
of equations to be solved by at least |(iV — n){N — n — 1) to give a full family of solutions. 
Thus we can expect solutions provided: 

d (nr{N - n) + i(iV - n){N - n + 1)) > i(iV - n){N - n - 1) 

+{d - 1) {nr-i{N -n) + |n^(n^ - 1) + nr{N - n) + |(iV - n)(iV - n - 1)) (42) 



Simplifying this gives 



N -n> 



nririr — 1) 



(43) 



2 - ^.-1) ' 

where we have made the assumption that the denominator on the right hand side is positive 
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Note the right hand side in (43) does not depend on the total dimension, n, of the 
matrices Ai, but just on rir-i and rir- Thus the size of the extensions needed in this special 
case may well be substantially smaller than in general. The application to cubature of this 
special form of commuting extensions will become clear in section 5.3. 



5 Cubature Formulae and Commuting Extensions 

Much of the contents of this section, with some additions, are also discussed in [12]. In 
section 5.1 we give a novel presentation on the subject of Gaussian quadrature; in section 
5.2 this is extended to the case of multidimensional cubature. Section 5.3 discusses some 
practical aspects and consequences of the results of 5.2. 



5.1 A Novel Approach to Gaussian Quadrature 

In Gaussian quadrature we wish to find q + 1 nodes Xq, . . . ,Xg and q + 1 weights Wq, . . . ,Wq 
such that the quadrature rule 



/ 



w{x)j\x)dx K,^Wif{xi) (44) 

' i=0 

is exact whenever f{x) is a polynomial of degree at most 2g + 1. Here is some interval 
or union of intervals and w{x) > is a suitable weight function. Throughout this paper we 
only consider quadrature and cubature rules with positive weights, i.e. iWj > 0. 

Note: In this subsection there is no mention of commuting extensions so we allow ourselves 
to break our index conventions. In this subsection alone indices i, j, /c run from to q. 

Denote by Vq the space of polynomials of degree at most q with the inner product 

{a\b) = w{x)a{x)b{x)dx , Ma.beVq . (45) 

Let Ilg be the projector from Vq+i onto Vq parallel to its orthogonal complement i.e. the 
obvious orthogonal projection onto Vq with respect to the inner product above. We define 
the operator X '■ 'Pq ^ 'Pq XP = ^qXp for all p E Vq. Since x is self adjoint there is an 
orthonormal basis of Vq consisting of cigcnfunctions {ui} of = Aj?7,j. Associated with 

X there is a symmetric bilinear form X on Vq defined by X{a, b) = (a|x6), or equivalently 

X(a, 6) = j w{x)a{x)xb{x) dx , Va, 6 G Vq . (46) 

X is diagonalised in the basis {ui\, X{ui,Uj) — A^Sij. 
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We prove below that the eigenvalues {Aj} of x provide nodes for a Gaussian quadrature 
formula of degree 2q + 1. Our treatment is the reverse of the classical presentation of 
Gaussian quadrature [1], [3], see the explanation after theorem 8. 

We need the following remarkable lemma: 

Lemma 2. (The S lemma) Let p be an arbitrary polynomial in Vq+i- Then 

{p\ui) = {l\uMK) , (47) 
i.e. the inner product of p with Ui is determined, up to normalization, by evaluation of p at 

Note: In this lemma p is allowed to be in Vq+i. 
ProoL We prove recursively for j that 

= (l|M,)Ai , j = 0,...,q + l. (48) 

For j — the statement is trivial. For j > 0, 

{x^Ui) = {xx^~^\ui) = {x^-^\xui) = Ai{x^-^\ui) . (49) 

This provides the recursive step proving (48), the full result follows by hnearity • 

X is an approximation of the operator x and this lemma shows that {wj}, the eigenfunctions 
of X, share with 6 functions, which are "eigenfunctions" of x, the property that projection of 
a function on either is done by its evaluation at the appropriate eigenvalue. This similarity 
is the reason for the name we give to the S lemma. 

With this lemma it is almost immediate to prove the main result we want: 

Theorem 7. Let / be a polynomial of degree at most 2q + l. Then 

/ w{x)f{x)dx = X^(l|ii,)V(AO , (50) 

i.e. the quadrature rule 

/ w{x)f{x)dx^j2i^\uiyf{K) (51) 

is exact of degree 2q + 1. 

Proof. Again we prove the result for /(x) — x^ , j — 0, . . . ,2q + 1, the full result follows by 
linearity. For j > 1, choose integers rii, ^2 between and q such that j — rii + n2 + 1. We 
then have 

/ w{x)x^dx = X(x"Sx"^) = X (^(x">fe)«fe,^(x"^|M,)K, ) = f](l|Ki)'A] . (52) 

''^ \k=0 1=0 / i=0 
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J2{l\ui)uA = Y.{l\ui)\ (53) 

i=0 / i=0 



In the last step we have used the 5 lemma twice. 
For j — observe that 

/ w{x)dx ^ ^ (^{l\uk)uk 

\k=0 

by the orthonormality of the Mj. • 

Theorem 7 relates the spectrum of x to the nodes of Gaussian quadrature. It is in fact 
easy to prove other facts in the theory of Gaussian quadrature using our approach. For 
example, to see that none of the weights vanish just take p — Ui in the 5 lemma. To see 
that when Q is a single interval [a, b] the nodes must be in its interior of the interval just 
observe that 

rb rh 

b — Ai = / w{x){b — x)ui{x)'^ dx > 0, a — Ai = / w{x){a — x)ui{x)'^ dx < . (54) 

Ja J a 

We also obtain the following widely known characterization of the nodes: 

Theorem 8. The nodes Aj are roots of any nontrivial degree g + 1 polynomial orthogonal 

to Vq. 

Proof. Let p be a nontrivial degree g + 1 polynomial orthogonal to Vq. Then using the 5 
lemma, = {p\ui) = {l\ui)p{Ai). Since is nonzero, this gives p(Aj) = 0. • 

Our presentation on Gaussian quadrature is the reverse of that in [1] , [3] . The starting 
point in [1], [3] is that the nodes in degree 2g + 1 Gaussian quadrature are roots of the 
degree g + 1 polynomial p from theorem 8; it is then shown that the eigenvalues of a matrix 
representation of x are equal to these. Here we have gone in the other direction; without a 
priori assumption of the existence of a Gaussian quadrature formula, we have shown that 
the eigenvalues of x s-^® quadrature nodes and as a consequence of the 5 lemma we also 
obtain the fact that they are roots of the degree g + 1 polynomial p from theorem 8. 

As we shall see in the next subsection, our approach to Gaussian quadrature allows a 
generalization to higher dimensions. However, in the generalization of the 5 lemma the 
polynomial p is restricted to Vq^ not "P^+i, and this means that while theorem 7 can be 
generalized, theorem 8 cannot, at least not immediately. So the characterization of cubature 
nodes as roots of a polynomial (or a set of polynomials) seems to be lost. However the 
characterization of nodes as eigenvalues persists, as we now set out to show. 

5.2 Generalization to Cubature 

(d + q\ 

We denote by Vq the n = I ^ 1 dimensional vector space of polynomials in d variables 
Xi, . . . ,Xd ol total degree up to q (total degree is defined by degree(x^^x^^ . . . x^'^) = mi + 
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1712 + .. ■ + rnd). 

An A^-point d-dimensional cubature formula 



N 

w{x)f{x)d'^X Waf{Xa) (55) 



is said to be of degree D if it is exact whenever f{x) is a polynomial of total degree at most 
D and non-exact for at least one polynomial of total degree D + 1. Here f2 is a suitable 
region in R'^ and w{x) > a suitable weight function. The weights Wa are assumed positive. 

We supply Vq with the inner product 

(alb) = I w{x)a{x)b{x) d^x , V a, 6 e P„ , (56) 
Jn 

and define 11^, the orthogonal projection operator from P^+i onto Vq with respect to the 
above inner product, in the obvious way. We can then define d self adjoint operators 
Xi,...,Xd on Vq by 

XiP = ^qXiP , \/peVq , (57) 
with related symmetric bilinear forms : Vq x Vq ^ R, 

Xi(a,b) = {a\xib) = j w{x)a{x)xib{x) d'^x \/a,b e Vq . (58) 

Generally \xi, Xj] 7^ so we can not find a basis of Vq in which all the Xi simultaneously 
diagonalised, and we do not have a direct analog of the one-dimensional case in which the 
eigenvalues of the single operator x served as quadrature nodes. We shall show, however, 
that there is a correspondence between cubature rules and spectra of certain commuting 
extensions of matrix representations of the operators xi, . . . , Xd- As a first step towards 
this we prove the following: 

Theorem 9. Let the n x n matrices Ai, . . . ,Ad be the representations of the operators 
Xi,---,Xd in an arbitrary orthonormal basis {Cq} of Vq (so {Ai)ab = (ealXi^b))- Suppose 
that for the region Q and weight function w{x) we have a degree 2g + 1, point cubature 
rule with positive weights. Then there exist N x N symmetric commuting extensions of 
Ai,...,Aa. 

Proof. Suppose the cubature rule takes the form 

N 

w{x)f{x)d'^X ^ ^af{Xa) ■ (59) 



/ 

JQ 



a=l 

Then, since all integrands are of degree at most 2g + 1, we have 

N 



Sab = {ea\eb) = w{x)ea{x)eb{x)d'^X = y2'>^aea{Xa)eb{Xa) , (60) 

N 

{Ai)ab {calxieb) = w{x)ea{x)xieb{x)d'^x ^y2'^c,ea{xa){xa)ieb{xa) . (61) 
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Define the n x N matrix Q by 




(62) 



and N X N diagonal matrices A 



1 



. , Ad with diagonal entries (Aj) 



aa 




(63) 



Equations (60)- (61) read 



nxn 



A, = QKQ^ . 



(64) 



Using theorem 5 in section 2 we conclude that Ai, . . . , have NxN symmetric commuting 
extensions. • 

It is natural to ask whether the matrix commuting extensions of theorem 9 arc repre- 
sentations of commuting extensions of the operators {xi} in some N dimensional space of 
functions V which includes Vq as a subspace. In other words, is it possible to extend the basis 
{Ca} of Vq to an orthonormal basis of V by adding N — n orthonormal functions e„+i, . . . , Cat 
in such a way that the NxN matrices (Aj)^^^ = {calxie/s) = J^w{x)ef^{x)xiei3{x) dfix com- 
mute? Unfortunately, in all but the simplest cases, we could not find, nor prove the existence 
of, such functions e„+i, . . . ,6^. However, even though we can not view the matrix commut- 
ing extensions of theorem 9 as representations of operator extensions of the Xj, they do 
satisfy a certain compatibility condition with the x« which we prove in theorem 10. 

Let us introduce the N dimensional vector space y(= R^), with the standard inner 
product, whose elements we denote in bold face. Vq is mapped to a subspace of V by the 
inclusion operator l : Vq ^ V which is defined by ie\ — ei,...,te„ = en, where {Ca} 
are the first n members of the standard basis of V . Extend {ea} to an orthonormal basis 
{go} ofV by adding any orthonormal basis {en+i, . . . , bn} of the orthogonal complement 
of span({ea}) in V . Even though our attempts to extend Vq with functions failed, now 
we are extending with A^-tuples. Define the obvious projection operator tt : V ^ Vq hj 
Trei = ei, . . . , yrSn = e„, TrSn+i = . . . = ttcn = E Vq] clearly ttl = I, the identity operator 

on Vq. 

Recall that the essential step in the proof of theorem 5 is extension oi Q to a,n N x N 
orthonormal matrix Q by appending any N — n orthonormal rows. In this way Ai, . . . , A^, 
NxN symmetric commuting extensions of the A^ are constructed in theorem 9, where 
Ai — QAiQ^. Since the Ai, . . . ,Ad mutually commute and are symmetric, there exist N 
orthonormal common eigenvectors Ua & V, such that AiUa — AiaUa- The matrices Ai are 
given in the basis {Bq,} and Q is the transformation between this basis and the eigenvector 
basis {Uq}. Note that the rows of Q give the coordinates of the vectors {Bq,} in the basis 
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{Uq,}, hence the extension of Q to Q by adding arbitrary N — n orthonormal rows is nothing 
but the extension of {ea} to {co,} the orthonormal basis of V described above. The reader 
is reminded that at present we are assuming the existence of a cubature formula hence the 
eigenvalues Aj^ are defined in (63). We can now state: 

Theorem 10. The commuting extensions of theorem 9, {^i}, satisfy the following com- 
patibility condition with the operators {xi}, 

Aitp = iXiP = iXiP, Vp e Vq-i , (65) 

where Vq-i is regarded in the natural way as a subspace of Vq. 

Note: Applying tt to (65) gives TrAjtp = XiPi P ^ 'Pq-i-i which is automatic as Ai is an 
extension of Ai. However, (65) contains more information than this, and is not true for an 
arbitrary extension of A^. 

Proof. We first prove that for any p &Vq and any eigenvector Uq, 

{ip\\lo) = Vw^p{Xa) , (66) 

where Wa, Xa, are the weights and nodes of the cubature formula whose existence is assumed 
in theorem 9. We already noted that the rows of Q from the proof of theorem 9 give the 
coordinates of the N vectors ea in the basis {uq}, in particular the rows of Q give the 
coordinates of ea = LCa, a — 1, ... ,n. Recall that Qaa — \^w^Gaixa), thus (66) is proven for 
basis elements Ca E Vq. The proof for general p E Vq follows by linearity. 

We now expand the left hand side of (65) in the basis {Ua}- For any p e Vq-i 

{AiLp\Ua) = Aia{ip\na) = KaV^Pi^a) = V^{^a)iP{Xa) , (67) 

using (66) and (63) respectively in the last two steps. To expand the right hand side of (65) 
note that XiP — XiP £ Vq for all p e Vq-i. Invoking (66) again we obtain 

{LXip\Ua) = \/w^iXa)iPiXa) , (68) 

which completes the proof. • 

Note that taking p = 1 in (66) gives Wa = (^IIuq)^. In theorem 9 wc saw that the eigenvalues 
of the commuting extensions are related to cubature nodes, here we obtain a relation between 
the weights and common eigenvectors. 

We shall see in section 5.3 that with an appropriate choice of basis the compatibility 
condition of theorem 10 implies that the commuting extensions have certain off-diagonal 
zero blocks; in particular this special structure aids computation of commuting extensions. 
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The obvious question to ask at this stage is whether there is a converse of theorems 9 
and 10. That is, suppose we have Ai, . . . ,Aci, N x N symmetric commuting extensions of 
Ai, . . . , Ad, which satisfy the compatibihty condition AiLp — iXiP — tXip, Vp e Vq-i . Can 
we build a cubature rule, without a priori assumption of its existence, using the eigenvalues 
and eigenvectors oi Ai, . . . , Ad7 In theorem 11 we give an affirmative answer to this. The 
treatment follows the presentation on Gaussian quadrature from section 5.1; in particular 
we start with a 5 lemma. Note that given the commuting matrices Ai we can find their 
diagonal representations Aj, but we do not assume in advance any connection of the Aj with 
cubature nodes. 

Lemma 3 (The multidimensional S lemma) Suppose the commuting extensions satisfy the 
compatibility condition A^Lp — iXip — iXiP for all i = 1, . . . , d, and for all p e Vq-i. Then, 
for any p EVq 

{ip\Ua) = p{K) , (69) 

where the points Aq G R*^ have entries {h-ia, ■ ■ ■ , A^q), all eigenvalues of Ai, . . . , Ad, satisfying 

^jUa — Ajo,1J.o;. 

Proof. We prove (69) for monomials For p = \ the statement is trivial. 

For any other monomial p &Vq we can write p = XiP, where p e Vq-i- Then 

{ip\Ua) = {iXip\Ua) = {Aiip\Ua) = {ip\AiUa) = Aia (i-plu^) . (70) 

Here the compatibility condition was used in the second step. Repeated application of (70) 
completes the proof for monomial p; the full result follows by linearity. • 

Note: Recall that we do not know how to relate y to a space of polynomials (or other 
functions) in a way which gives commuting extensions of the operators Xi, ■ ■ ■ , Xd- In par- 
ticular, we can not interpret the eigenvectors u^^ as polynomials (or other functions). Thus 
our present state of understanding allows us to view the Uq, as "5 vectors" in V and not as 5 
functions, which was possible in the 1-dimensional case. Moreover we can not identify Vq+i 
with a subspace of V. Hence, in contrast to the one dimensional case, we restrict p G Vq in 
the multidimensional S lemma thereby losing the immediate connection between cubature 
nodes and roots of polynomials in Vg. 

We are now fully prepared for the converse statement to theorems 9 and 10: 

Theorem 11. Let Ai, . . . , Ad be the representation of the operators Xi, • • • , in an or- 
thonormal basis {ca} of Vq. Let Ai, . . . , Ad he N x N symmetric commuting extensions of 
Ai, . . . ,Ad satisfying the compatibility condition A^Lp = iXiP = iXiP for all p e Vq-\. Then 
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for every polynomial / in 7^2^+1 



N 



[ W{x)f{x)d'x = J2{a\Uaff{Xa) . (71) 

Here the Uq, are joint eigenvectors oi Ai, . . . ,Ad, satisfying A^Ua — AjaU^, and the points 
Xa e R*^ have entries (Ai^, . . . , Ada)- 

Proof. Recall the symmetric bilinear forms on Vq defined in (58). Given the commuting 
extensions, we can introduce symmetric bilinear forms Xi on V defined by Xi{u,v) — 
{u\AiY). Since the Ai are extensions of the Ai we have X{pi,p2) — Xi{Lpi, LP2) for all pi,p2 
in Vg. The Xi are simultaneously diagonalized in the basis {u^}, Xi{Ua, up) = AiaSap- 

It is sufficient to prove the statement of the theorem for monomials. For / = 1 we have 



/ w{x)d''x = (1|1) = = /X^(tl|u„)u„ 

-'^ \a=l 



N \ N 

x:(^l|u^)u,) = ^(.l|u„)^ (72) 

13=1 I a=l 



by the orthonormality of the Uq,. Note that in the second expression in (72) the inner 
product is taken in Vq., in subsequent expressions it is taken in V . 

Any other monomial in V2q+i can be written in the form / = 2:1/1/2 for some monomials 
/i, /2 £ Vq and some i. Note that use of the multidimensional S lemma is possible since we 
assume the Ai satisfy the compatibility condition, and so 



f w{x)f{x)d'^X = [ w{x)fi{x)Xif2{x) d'^X = Xi{fij2) 



(TV TV 
^(t/l|Ua)Ua, J2{'^M'^f3)Uf3 
a=l /3=1 
N TV 
= ^ Ai„(t/l|Ua)(t/2|u«) = X!(''l|"«)^Ai„/l(A«)/2(Aa) 
a=l a=l 
N 

= i:(.l|u„)V(Aa) . (73) 

a=l 

Thus (71) is proven for monomial /; the full result follows by linearity. • 

Theorems 9,10,11 give the main result of this paper, that N point, odd order cubature 
formulae with positive weights are equivalent to symmetric commuting extensions, satisfying 
the compatibility condition, of matrix representations of the operators xi, . . . , Xd- 



5.3 Discussion and Consequences 

Our findings give a new computational approach to the derivation of cubature formulae. If 
appropriate commuting extensions are numerically found their simultaneous diagonalisation 
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will give the cubature rule in (71). To numerically obtain the matrices {Ai} we introduce 
an orthonormal basis of Vq consisting of an orthonormal basis of Vq (a constant function ei 
with ||ei|| = y^(ei|ei) = 1), extended to one of Vi, extended to one of V2 etc., i.e. a basis 

fd + q\ 

{ca}, a = 1, . . . ,n = [ ^ 1, such that 

Ci is an orthonormal basis of Vq 

ei, . . . , 6(1+1 is an orthonormal basis of Vi ^^^^ 

ei, . . . , ei((^^i)(^_,_2) is orthonormal basis of V2 

etc. 

Such a basis can be obtained from the monomials {x^^ . . . x^'^}, nii + . . . + ma < q, hy the 
Gram-Schmidt procedure. Note that all basis elements Ca of degree m or more are orthogonal 
to Vm-i- This choice of basis and the fact that {Ai} represent the operators {xi} imply 
that the A^ have tridiagonal block form as in (36), with q + 1 blocks on the diagonal. In 

the notation of section 4, rii = 1 and for m = 2, 3, . . . , g + 1, = dimP^.i — dim'Pm-2 

(ni = 1, n2 = d, ris = d{d + l)/2 etc). Moreover, the commutator of any pair of Aj's 

/q + d-l\ fq + d-l\ 
is zero apart from a single block of size x in the bottom right 

V d-1 J \ d-1 J 

fq + d-l\ 

hand corner (note = rig+i = dimP^ — dimVq-i). The compatibility condition, 

V d-l J 

together with the fact that en+i, . . . , Bn are orthogonal to iVq, imply that the bottom left 
(N — n) X dim{Vq-i) block of Ai is zero and by symmetry so is the corresponding block in the 
upper right. Thus the commuting extensions {Ai} we seek are precisely those in tridiagonal 
block form as in section 4. Note also that since our first basis element ei is a constant 
polynomial the cubature weights are obtained from the first entries of the eigenvectors Uq,, 

In the case d — 2, we note that the matrices Ai, A2 have n^-i = q and rir — q+1 in the 
notation of section 4, and thus from (43), which is based on counting degrees of freedom, 
the expected size of the commuting extensions is 

N>n+'-^. (75) 

Using n — dimP^ = |(g + l)(g + 2) we obtain 

(2q + 2)i2q + 3l 1 
N > = -dimP2g+i ■ (76) 

This is exactly the number of nodes we expect from counting degrees of freedom in a 2- 
dimensional cubature formula of degree 2q + 1 . 
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For d > 2 the Aj_ have a more subtle structure that requires refinement of the discussion 
leading to equation (43). However, the equivalence of cubature formulae and commuting 
extensions (satisfying the compatibility condition) allows us to estimate N by easily counting 
degrees of freedom in a general d-dimensional cubature formula of degree 2q-\-l. Thus, 

N > ^dimP2,+i . (77) 

We emphasize again that such calculations are not rigorous and the inequalities obtained in 
this way can serve only as recomendations for choice of N, indeed certain cubature formulas 
with less points are known. 

In section 6 we shall give some first examples of computation of cubature nodes using our 
approach. But before we do this we present two theoretical consequences of the equivalence 
between cubature formulae and commuting extensions. 

Theorem 12. Let N be the number of nodes in a degree 2g + 1, d-dimensional positive 
weight cubature rule. Then 

fd + q\ 1 

A^>f ^ 1 +-maxi,,rank([A,A,]) , (78) 

where Ai, . . . ,Ad are the matrix representations of the operators xi, . . . , Xd on Vq. 

Proof. By theorem 9 an N point cubature rule gives N x N commuting extensions of 
the matrices Ai. By theorem 2, section 2, the size of such extensions is at least dimP^ + 
|maxj jrank( [A^ ,Aj]). • 

Notes: (1) As mentioned in the introduction, theorem 12 has its origins in the work of Moller 
[6]. A statement of the result in a form that clearly corresponds to our statement can be 
found in [17], which cites [18] and [19]. Our proof, however, is a substantial simplification. 
(2) It is informative to compare the lower bound of theorem 12 with estimates based on 
parameter counting. As a consequence of our previous remarks on the block structure of 

so the second term in (78) is typically a small fraction of the first term. Consequently for 
large q the right hand side of (78) is much smaller than the number of nodes we expect from 
counting degrees of freedom in a degree 2g + 1, d-dimensional cubature formula, which is 

This comparison indicates why the lower bound on the number of points needed for a 
cubature formula is rarely attained. 
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Theorem 13. Let Ai, . . . ,Adhe matrix representations oi Xi, ■ ■ ■ ,Xd- In any degree 2g + 1, 
d-dimensional, positive weight cubature rule, and for each i, there is a node Xa with {xa)i 
less than or equal to the smallest eigenvalue of Ai, and a node xp with (x/j)^ greater than 
or equal to the largest eigenvalue of Aj,. 

Proof. By theorem 9 a cubature rule of degree 2q + 1 gives commuting extensions of 
the matrices Ai with the nodes composed of the eigenvalues of the extended matrices. By 
theorem 6 in section 2 the smallest /largest eigenvalue of the extended matrices is less/greater 
than or equal to the smallest /largest eigenvalue of the matrices before extension. • 

Note: As far as we are aware this theorem is not even known ior d — 1. For d — 1 the 
theorem says that any A'"-point, positive weight, degree 2q + 1 quadrature rule must have 
a node less/greater than or equal to the smallest /largest Gaussian quadrature node. Thus 
Gaussian quadrature has the property that the span of the nodes is the smallest possible, 
amongst all positive weight quadrature rules, with any number of points, that are exact to 
the same degree. 



In this section we briefly discuss 2 examples of computing cubature rules via commuting 
extensions, both are in 2 dimensions. In section 6.1 we consider the classic question, first 
studied by Radon [20], of finding 7 point cubature rules which are exact for polynomials 
up to degree 5 (dimT's = 21 for d = 2). This involves 7x7 extensions of a pair of 6 x 6 
matrices of the tridiagonal block form of section 4, and we provide a reliable algorithm for 
this. In section 6.2 we present some new cubature rules for integration on the entire plane 
with weight function w{x,y) = e~^^~^^. These were computed by commuting extension 
techniques, though, as explained in section 3, our current algorithms for this are poor, so 
we give limited details. 

6.1 Radon type formulae 

Given a region on the plane and a suitable weight function w{xi, X2), a Radon formula is 
a cubature rule of the form 



6 Examples 




a=l 



(81) 



which is exact for all polynomials of degree no more than 5. 
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To construct such formulae we proceed as follows: First choose a basis ei, 62, 63, 64, 65, eg 
of V2 which is orthonormal with respect to the inner product 



{a\h) 



a, 



w{xi, X2)a{xi, X2)b{xi, X2) d?x , 



(82) 



and such that ei is of degree 0, 62, 63 are of degree 1, and 64, 65, Cq are of degree 2. Typically 
we do this by applying Gram-Schmidt orthonormalization to the basis 1, xi, X2, xf, X1X2, x\. 

Having constructed the orthonormal basis we compute the matrices A\ and A2 via: 



{A 



i)ij 



w{xi, X2)ei{xi, X2)xiej{xi, X2) d?x , 



X . 



(83) 
(84) 



By orthogonality we will have (Ai)ij = (Ai)ji = (^2)11 = (^2)41 = for i = 4, 5, 6 (this is 
the meaning of "tridiagonal block form" in this case), and the commutator C = [Ai, A2] will 
be all zero except for a single antisymmetric 3x3 block in the lower right hand corner. The 
commuting extensions Ai, A2 should take the form given in (18), where a, b are 6-dimensional 
column vectors, with the first 3 entries vanishing (these arc the dimP^_i x [N — n) zero 
blocks described in section 5.3), and a,j3 scalars. Following the arguments in lemma 1 from 
section 2, a,b,a,j3 must satisfy (19)-(20). To solve these equations wc proceed exactly as 
in the lemma. First we construct a specific pair of 6-dimensional vectors v, w satisfying 



[Ai, A2] + vvF - wv^ = 



(85) 



Because of the special structure of [Ai, A2], v,w can be taken to have zeros in their first 
3 entries, and determining the other entries is equivalent to finding 2 vectors in with a 
given cross product. Clearly v, w are linearly independent and give a basis for Im([74i, A-^). 
Once such v,w have been found, (19) gives 

a = \v + iiw , b = pv + pw , (86) 

where Xp — fiu — 1, and (20) gives the requirement 

(/9A — ai')v + {(3/1 — ap)w + vA\v + pAiw — \A2V — pA2W — , (87) 



or 



V w AlV Aiw A2V A2W) 



/ f3X- ai'\ 
(5p — ap 

V 

P 

-A 
\ -p ) 



(88) 
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Since v,w have zeros in their first 3 entries, and (Ai)ij = (^2)11 = for i = 4,5,6, the 
entire first row of the matrix {v w Aiv Aiw A2V A2W ) is zero. So it is singular and 
a nontrivial solution oi {v w Aiv Aiw A2V A2W ) A; = is guaranteed. To complete 
construction of a commuting extension all we need is to find a, P, A, /i, u, p with Xp — /ii/ — 1, 
and such that the column vector in (88) is in the kernel of the matrix in (88). Using 
the freedom to rescale vectors in the kernel, it is straightforward to conclude that given 
a nontrivial vector k in the kernel of the matrix, there will be an associated commuting 
extension if and only if k^ke — k^k^ > 0, and we must take 

A = — c/cs , p, — —ck^ , u — cks , p — ck^ , 

a = c^{k2k5 - kike) , P = c^{kik4 - k2kz) , (89) 
where c = , ^ , , 

V«3«6-«4«B 

To summarize: to find a Radon formula one should: (1) Construct the orthonormal basis 
{ca}. (2) Construct the matrices Ai,A2. (3) Find 6-dimensional vectors v,w with top 3 
entries zero, and such that [Ai, A2] + vw'^ — wv'^ = 0. (4) Compute the kernel of the matrix 
( V w AiV AiW A2V A2W). Then for each vector k in the kernel with ksk^ — k^k^ > 
there is an associated commuting extension given by (18), where a,b are given by (86) and 
a, (3, A, p, u, p by (89). The commuting extensions can be simultaneously diagonalized, using 
the algorithm in [15], to obtain the nodes and weights of the cubature rule. 

This procedure should be compared with Radon's original work [20]. In practice we 
have not found a case in which there is a vector in the kernel with k^k^ — k^k^ < 0, but 
have no explanation why this is so. Generically the kernel is one dimensional and there is a 
unique Radon formula. For the case of equal to the circle or the square each with uniform 
weight w{xi,X2) — 1, the kernel has dimension 2, giving a one parameter family of Radon 
formulas. In the case of the square with vertices (—1, —1), (—1, 1), (1, 1), (1, —1), all the 
Radon formulae have a single node at the origin with weight | and 3 pairs of diametrically 
opposed nodes on the circle xl + x^ = j^. (In the literature there is often only mention of the 
case of Radon formulae for the square when one pair of nodes lies on one of the coordinate 
axes.) In the case of the unit circle, there is a single node at the origin with weight | and 
3 pairs of diametrically opposed nodes on the circle x^ + X2 = ^ ; in this case there is full 
rotational symmetry, and the different formulae are related by rotation. 

As an explicit example of nonstandard Radon formulae, we give here some results for 
the square with vertices (—1, —1), (—1, 1), (1, 1), (1, —1), with a square with vertices (| — 
r, I - r), (I - r, I + r), (| + r, | + r), (| + r, | - r) removed. Here < r < | and the 
weight is uniform. Figure 1 displays the location of the nodes for r = ^, i = 0, 1, . . . , 8, 
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along with the boundaries of the squares removed, and the circle + x| = ^ on which 
the nodes lie for r = 0. For r = we in fact plot the nodes for a very small nonzero value 
of r, so as to have a unique result. For values of r up to | all the nodes are inside the 
domain of integration. For r — ^ and r — one node lies inside the small square that is 
removed to obtain the cubature domain. For r = | one node lies outside the larger square, 
at {xi,X2) ~ (0.1844, 1.0360), but has low weight (about 3.25% of the sum of the weights). 

Recall that our approach for constructing cubature formulae involves first constructing 
commuting extensions, and then diagonalizing them. For the rational values of r discussed 
in the previous paragraph, the commuting extensions can be constructed exactly, and we 
do not need to apply the numerical algorithms for commuting extensions from section 3 
(the S{Q, Aj) algorithm does not work particularly well for the matrices in this subsection). 
Even though the commuting extensions can be found exactly, finding the cubature nodes 
(the eigenvalues of the extended matrices) can only be done numerically. 

6.2 Gaussian weight on the plane 

As already explained, to construct degree 2g+l cubature formulae on a region Q in the plane, 
we need to find symmetric commuting extensions (satisfying the necessary compatability 
conditions) of the nxn matrices defined in (83)- (84), where n — dimT^^ = (g + 1) Qgr + 1^. 
By counting degrees of freedom, we expect A'"-dimensional symmetric extensions if > 
(q + l) + l), see the discussion leading to (76). We work in a basis of the type discussed 
in subsection 5.3, so the matrices are all of the tridiagonal block form of section 4. Using the 
notation of section 4, we extend each matrix by adding 2 blocks, of size (q + l) x (N — n), 
and CKj (symmetric), of size {N — n) x {N — n), here i — 1, 2. Thanks to the freedom in choice 
of extensions described in theorem 3, one of the symmetric added blocks ctj can be chosen 
diagonal. Taking this into account we expect N — n > ^q{q + 1). Even so, the number of 
entries added in the blocks Oi, 02, cti, 02 grows as g^. 

For applications in quantum mechanics, cubature on the plane fl = with weight 
function w{xi,X2) = e~^i~^2 is important. For this case we have computed degree 11, 
13, 15 and 17 cubature formulas with 26, 35, 46 and 57 nodes respectively, by simultaneous 
diagonalization of commuting extensions. These were found using the gradient fiow approach 
mentioned at the begining of section 3. For degrees 11 and 13 we have also succeeded to 
compute the necessary commuting extensions using the S{Q,Ai) algorithm on which section 
3 focused. To take into account the zero blocks of the commuting extensions an extra term 

E E E iiQKQ^)a^) (90) 

i=l a=l b=n+l 
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-0.6 



-0.8 



Figure 1: Radon formulae on a square with a smaller square removed, with uniform weight. 
As the size of the smaller square tends to zero, one node tends to the origin and 6 to the 
circle x\ + = As the size of the smaller square is increased, one node moves outside 
the domain of integration, and eventually out of the larger square. 
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Figure 2: Nodes for a 46 point, degree 15 (left) and a 57 point, degree 17 (right) cubature 
formula for the plane Q = with Gaussian weight w{xi, X2) = e~^i~^'2 

was added to S of (32), this gave a substantial improvement in performance (compared 
to the same algorithm with S not including the extra term). The nodes and weights are 
available at 

http : //www. math. biu. ac . il/~schif f /commext .html . 
Figure 2 displays the location of the nodes in the degree 15 and 17 formulae. The formulae 
found are not necessarilly minimal, and are certainly not unique, because of rotational 
symmetry. It seems our degree 13, 15 and 17 formulae are new, the one with order 17 
having a smaller number of nodes compared to existing formulae of the same order, see [21]. 

7 Summary and open questions 

The central results of this paper are theorems 9,10,11, which prove the equivalence of cuba- 
ture formulae and commuting extensions satisfying the compatibility condition (equivalent 
in an appropriate basis to requiring certain zero blocks in the extension matrices). This 
raises the questions of existence and methods of computation for commuting extensions. 
Our knowledge of the theory of commuting extensions is summarized by theorems 1 to 6, 
and in section 3 we have described our initial attempts at their computation. 

There is clearly enormous potential for further work here. In the context of our main 
topic, the connection between cubature formulae and commuting extensions, there is one 
nagging question that we have indicated several times in section 5: The vector space V 
on which the commuting extensions act does not yet have an interpretation as a space of 
functions (or maybe even polynomials). For numerical work in quantum mechanics it would 
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be a major advantage if we could contstruct finite dimensional function spaces containing 
the space of all polynomials of a given degree as a subspace, on which the natural projections 
of the operators Xi commute. The existence (or nonexistence) of such spaces is a topic we 
hope to investigate, see also [5] . 

Another question left open in our work is that we have not given an existence proof 
of cubature formulae from the commuting extension viewpoint. Although theorem 1 in 
section 2 guarantees the existence of commuting extensions of an arbitrary set of matrices, 
it does not guarantee extensions in the form we need to be able to apply theorem 11. An 
existence proof for commuting extensions of the required form would provide an alternative 
approach to Tchakalojf 's theorem [22] that guarantees (for any suitable domain fl and weight 
function w{x)) the existence of positive weight cubature rules that are exact for certain sets 
of functions. 

The numerical question of computing cubature formulae is now subsumed under the more 
general question of computing commuting extensions, likewise the open sore that there is no 
good way to predict the minimal number of points needed for a cubature rule is subsumed 
under the question of finding the minimal dimension for commuting extensions. There are 
a number of points in the theory of commuting extensions which we feel may be improved, 
for example, existence of symmetric commuting extensions for symmetric matrices may 
well be provable, but we suspect the question of minimal dimension is extremely difficult. 
Fortunately, just because it is difficult theoretically does not mean answers cannot be found 
numerically, and we are hopeful that good algorithms can be devised that find commuting 
extensions of a given dimension, if they exist. The determining equations are linear and 
quadratic, and although there surely will be ill-conditioning in certain cases, it is hard to 
see why this should be so in general. We suspect that the poor performance of the S{Q, Aj) 
algorithm in section 3 is to do with the fact we inforced equation (34), and only searched 
over Q. Likewise gradient and Newton algorithms possibly give exaggerated importance to 
linear components of the system. There is much more that can be tried here. 

Another aspect to be considered in construction of cubature formulae/commuting exten- 
sions is symmetry in the domain D, and weight w{x). This will clearly influence the matrices 
Ai, which represent the natural projections of the operators Xi, and should be respected in 
construction of the extensions Ai, see also [5]. 

We hope very much that more applications will emerge for the notion of commuting 
extensions. The idea that noncommutativity can be resolved by introducing extra dimen- 
sions is a very natural one. In fact, we suspect that, more than the ranks or the norms of 
commutators, the size of minimal commuting extensions is probably the best measure of 
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how noncommuting a set of matrices is. 

The minimal size issue appears in other settings too. For example, given a set of m x n 
matrices Ai we can ask what is the smallest N such that there exist anmx N matrix U and 
a,n nxN matrix V, both with orthonormal rows, and N xN diagonal matrices A^, such that 
Ai = UAiV'-'-'. In our context, this provides a natural generalization of the singular value 

decomposition of a single matrix, in the same way that (26) provides the generalization of 
diagonalization of a single symmetric matrix. 
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